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Abstract. We consider the problem of imposing radiation conditions at artificial boundaries for the 
numerical simulation of wave propagation. Our emphasis is on the behavior and analysis of the error which 
results from the restriction of the domain. We give a brief outline of the theory of error estimation for 
boundary conditions. Use is made of the asymptotic analysis of propagating wave groups to derive and 
analyze boundary operators. For dissipative problems this leads to local, accurate conditions, but falls 
short in the hyperbolic case. A numerical experiment on the solution of the wave equation with cylindrical 
symmetry is described. We give a unified presentation of a number of conditions which have been proposed 
in the literature and display the time dependence of the error which results from their use. The results are 
in qualitative agreement with theoretical considerations. We find, however, that for this model problem it 
is particularly difficult to force the error to decay rapidly in time. 
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1. INTRODUCTION 


1. Introduction. It is a common feature of many wave propagation problems that 
they are naturally posed on unbounded spatial domains and that disturbances are radiated 
to the far field. The numerical solution of such problems requires their restriction to a finite 
computational domain, usually via the introduction of an artificial boundary. The choice 
of radiation boundary conditions there can have a dominant influence on the solution, 
particularly for long time simulations. For this reason, the study of numerical radiation 
conditions is an important aspect of the numerical analysis of wave propagation, and has 
been the subject of a great deal of research in recent years. In this paper, we discuss some 
aspects of the problem with an emphasis on the behavior of the error which results from 
introducing the artificial boundary and on the use of asymptotic analysis in its estimation. 

The outline of a general theory of artificial boundary conditions is fairly straightforward 
to develop. Consider a linear problem on a cylindrical domain: 

(!) L(d t ,d x ,x)u = f, x = (xi,x 2 , . . .) € [0,oo) x Q , t > 0, 

supplemented by appropriate boundary and initial conditions. (Note that coordinate sys- 
tems can be chosen so that the far field of exterior domains is transformed into a cylinder.) 
Suppose an artificial boundary is introduced at x\ = r and the radiation condition: 

(2) B v — g, 

is imposed on the approximate, finite domain solution v. The error resulting from the 
restriction of the domain, e = u — v, satisfies: 

(3) Le = 0, x G [0, r] x Q, t > 0, 

(4) Be = Bu — g, x [ = r, 

in addition to homogeneous initial conditions and homogeneous conditions at <9f2 and at 
x i — 0. To estimate ||e|| we need appropriate generalizations of the usual concepts from 
numerical analysis, stability and consistency. Stability in this case is simply the well- 
posedness of (3- 4), which we take to be the existence of a finite K independent of the data 
Bu — g such that, 

( 5 ) || e ll < K\\Bu — </||. 

Here the particular norms are determined by the differential operators L and B. In order to 
take account of any decay of the error off the artificial boundary, use of weighted norms may 
be of interest. Stability can be established, in principle, by studying Kreiss type conditions. 
Given stability, we see that the error will be small if ||J3u - g|| is; that is if the solution of 
the original problem approximately satisfies the radiation condition. Carrying further the 
analogy with the theory of discretization error, we may consider a sequence of approximate 
problems depending on a parameter e — *■ 0 such that, 

( 6 ) ||S(e)u - 5 (0|| = O ((ff(0) _1 ), 

which is simply a consistency condition. Clearly, (5) and (6) together imply the convergence 
of || e j| to 0. 
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2. ASYMPTOTIC ANALYSIS AND GROUP VELOCITY 


The problem of designing effective approximate radiation conditions is, then, to find 
B and g such that the stability condition is satisfied and || Bu — <jr|| is small. There is the 
additional consideration that the finite domain problem for v be reasonably easy to solve, 
which rules out arbitrary nonlocality in space and time for the operator B. The error 
analysis outlined above, though not typically discussed in papers on numerical radiation 
conditions, has been carried out in a few cases. Halpern and Rauch [11] consider general 
conditions for the wave equation which are exact in the high frequency limit and show that 
the error is small if the average frequency present is high. Bayliss and Turkel [2] use a long 
range asymptotic expansion of solutions of the wave equation to derive error estimates in 
powers of r” 1 . In joint work with Hariharan [10] we generalize these to convective wave 
equations. In all the cases above, as will be seen in the numerical example later on, the 
error estimates are nonuniform in time. Error estimates in r 1 are also given in [9] for 
dissipative problems. These will be discussed in more detail in the next section. 

2. Asymptotic Analysis and Group Velocity. To make further progress we sup- 
pose that the coefficients of L are independent of xi and t for x\ > r. Then, using Laplace 
transforms and eigenfunction expansions, u may be represented as a superposition of normal 
modes, 


( 7 ) 


u 


Y2 jf c i{p)qi(t - p< x )dp, xi>r 


(8) «(t, «c) = f e at+X ^y t (x 2l . . . ; s)ds. 

v ’ 2m Jc 

(For further details of this representation and the following arguments see [9].) If L is second 
order in x\, an exact boundary condition satisfied by u is given by: 

dui „ 

(9) ^ 

The generalized dispersion relation, A;(s), which defines the exact boundary operator for 
the /th mode, is difficult to obtain in analytic form except for problems with constant 
coefficients. Furthermore, (9) defines an operator which is nonlocal in both space and time 
and, hence, which violates the condition that the reduced domain problem be easily solved. 
To get an operator which is local in time, we must approximate A/(s) by a rational function 
in s. (Typically, locality in space is also obtained by restriction to a small number of modes 
or approximating the / dependence.) Let the approximate condition be expressed by. 

( 10 ) — = 

Then the transform of Bu, which must be made small to get a good error estimate, is given 
by: 

(11) Bui = (A/(s) - 6/(s)) iii. 

One way to make (11) small would be to make A; - 6/ small. However, we cannot expect 
to be able to enforce this everywhere in s space, as a low degree rational approximation 
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2. ASYMPTOTIC ANALYSIS AND GROUP VELOCITY 


to an arbitrary function is, in general, only accurate in restricted regions. We must, then, 
choose regions in s space where A / is to be approximated. Clearly, these should correspond 
to regions where ti/ is large. It is in identifying these that asymptotic analysis can play an 
important role. We note that the complementary problem of determining the stability of 
rational approximations is taken up in detail for a particular case by Trefethen and Halpern 
[ 12 ]. 

Asymptotic expansions of u/ for r large can be computed by the method of steepest 
descent and lead to a study of wave groups. These are defined by: 

(12) *{(,)= - 7> 

where 7 is real and nonnegative. We suppose that ( 12 ) defines curves, s = s*( 7). To derive 
boundary conditions, we must further restrict attention to small subintervals of these curves. 
We note that each wave group has an exponential decay rate, /i/( 7), given by: 

(13) W(7) = -*(**(7)7+ A,(«*(7))). 

We expect the large behavior of the solution to be dominated by wave groups with 
minimum decay rates, if they exist. Differentiating /i; with respect to 7 yields: 

(14) n\{ 7 ) - -S( 0 . 

This leads to the consideration of two typical cases: 

Case 1: The curve s*( 7) coincides with part of the imaginary s axis. 

From ( 14 ) we see that the decay rates are independent of the group velocity. This indicates 
that all wave groups may be important in the far field and, therefore, provides no guidance 
to the local approximation of the dispersion relation. Case 1 holds for hyperbolic problems. 
Consider, for example, the dispersive wave equation, 

( 15 ) Uft — u xx + u = 0, 
which has the dispersion relation: 

( 16 ) A = -(1 + s 2 )%. 

Solving (12) we find that curves of real, positive group velocity are given by the portions 
of the imaginary axis above i and below where the decay rates are /i = 0, and the 
nonnegative real axis where /j is positive. It is worth noting that the question of finding 
accurate and efficient boundary conditions for this simple problem has never been satisfac- 
torily settled. It has recently been taken up again by Engquist and Halpern [4] who propose 
an approximation to A which is accurate for s large and at s = 0. 

Case 2 : The curve(s) $*(7) cross the imaginary s axis. 

In this case ( 14 ) implies a local minimum (or maximum) of the decay rate. An accurate 
asymptotic boundary operator can then be constructed by approximating the dispersion 
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2. ASYMPTOTIC ANALYSIS AND GROUP VELOCITY 


relation near the crossing point. In [9] we consider a linear approximation which leads 
to particularly simple boundary operators involving first order derivatives in space and 
time. Assuming a global minimum for the decay rate, we were able to show that the 
error is proportional to r” 1 . Applications of this technique to simulations of diffusion and 
advection by a parallel flow [7] and the propagation of disturbances to Poiseuille flow [8] 
have illustrated its effectiveness. Case 2 typically holds for problems with dissipation. The 
simplest example is the advection-diffusion equation in one dimension: 


(17) 


U t “h XI x W xx — 0 , 


with dispersion relation: 

(18) A = 1 (l - (1 + 4*)*) . 


The curve of real group velocity coincides with the real s axis to the right of s = — so 
that the decay rate has a minimum at s = 0. A boundary condition, 6, which agrees with 
A in a neighborhood of 0 will be accurate. 

It is of interest to examine the connection between the two cases by considering a dis- 
sipative problem which approaches the dispersive wave equation as a dissipation parameter 
approaches zero. We take, as a particular example, 


(19) 


d 

at 



d 2 u 2 

u - —z + ct l u - 0. 
ox 1 


The dispersion relation has two branches: 

(20) A± = --L- (l + 2 vs + 2 u 2 a 2 ± (1 + 4 t / s )*) 2 . 


For v small A + is large and is associated with rapidly decaying wave groups. It is wave 
groups defined by A_ which are most important. We find two curves of real group velocity. 
The first coincides with the real s axis to the right of — The minimal decay rate, A_(0), 
is —a. A second curve can be found numerically. Its ends are the branch points of A_, 
$± = — j/a 2 ± ia. Following the curve from s . (_, the real and imaginary parts of s increase 
as we cross the imaginary axis. Eventually, the imaginary part reaches a maximum and 
decreases until the curve crosses the real axis. Expanding A_ for v small we find that the 
crossings take place at s — ±yj^ia -h 0{v). A local approximation to A, is given by: 

(2i) A_ ~ " ( ± ^ + T w2 + ^ (st ' 


We see that the decay rates are 0(v) so that these waves decay far less rapidly than the 
branch through s = 0. To derive boundary conditions which match the behavior of A_ at 
both crossings we may either consider products of first order conditions [9] or direct rational 
approximations. Following the latter approach we find: 


( 22 ) 


b 


4 s 2 + 9i/ , c* 2 s + 3a 2 
\/l2 s 
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3. THE CYLINDRICAL WAVE EQUATION: A CASE STUDY 


In practice, modifications of this condition to remove the singularity at 0 may be necessary. 
The corresponding differential operator is: 


(23) 


„ ( r— d 2 d 2 7 d 2 

B = V 12—— + 4— + 9t/a 2 — + 3a 2 
1 dxdt dt 2 dt 


Of course we require two conditions for a fourth order equation, but the additional condition 
is easily found. (For example, by differentiating the first.) It is tempting to use the v = 0 
limit of (22) for the dispersive wave equation. However it must be noted that the asymptotic 
error estimates are not uniformly valid as v — ► 0. The curves of real group velocity steepen 
and approach tangency to the imaginary axis. This implies that the decay rates are nearly 
constant. Nonetheless, it would be interesting to experiment with these conditions for v 
small, and we plan to do so in the future. 


3. The Cylindrical Wave Equation: A Case Study. As we have seen, the general 
asymptotic analysis of propagating wave groups described above does not yield a technique 
for computing accurate, local radiation boundary conditions in the hyperbolic case. However 
a number of results have appeared in the literature related to the accuracy of radiation 
conditions for the wave equation. As mentioned earlier, the work of Bayliss and Turkel 
[2] for exterior domains contains error estimates in r -1 . These are nonuniform in time. 
More recently, Engquist and Halpern [4], [5], Barry, Bielak and MacCamy [1] and Bielak 
and MacCamy [3] have considered the long time behavior and introduced modifications of 
conditions based on high frequency asymptotics to ensure the decay or at least suppress the 
growth of the error. 

To better understand both the long and short time behavior of the error for various 
boundary conditions, we have considered the following simple test problem: 


(24) 


d 2 u 

Tt 2 


1 <9 du 

(r — ), r > 1, 

r\ V / f — I 

r or or 


t > 0 , 


(25) 


. , cos2?rf , du 

u(M) = 1 - — — T , u(r,°) = — (r,0) = 0. 


An exact condition satisfied by u is: 

du sKq(ts) 


(26) 


dr K 0 (r.) “ = 


Here 2 = rs. All boundary conditions we consider may be written as approximations, 6(z), 
to T . For r large we must consider the behavior of T for z large: 

(27) T = - Z z ► oo . 

The long time behavior, however, corresponds to the limit z — ► 0 where we have: 


(28) 



z — ► 0 . 
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3. THE CYLINDRICAL WAVE EQUATION: A CASE STUDY 


We now list the conditions we have used in the experiment and compute their corre- 
sponding 6. 


Condition 1: 
(29) 


dv dv 
dr dt 


= 0. 


This simple condition follows, for example, from adding a correction to condition (32) to 
make 6(0) = 0 as suggested in [4]. The properties of 6 are: 


(30) 


b(z) = - z , 


(31) 


b-T = 




Condition 2: 

(32) 


dv dv v 
+ ~dt + 2r 


This condition follows from the high frequency analysis of Engquist and Majda [6] and is 
the first in the hierarchy of conditions proposed by Bayliss and Turkel [2], The properties 
of its 6 are: 

(33) b(z) = -z-^, 


(34) 


b-T = 


0(z *), z —> oo 
-i z 0 


Condition 3: 

( d d 5 \ / 9 d 1 \ 

(35) \T,' ¥ m + Tr)\Vr + ^i*Y r ) v - a 


This is the second condition in the Bayliss-Turkel hierarchy. It is chosen to match more 
terms in the large z expansion of T . The results of [2] imply, for fixed time, an error 
proportional to r -1 for Condition 2 and proportional to r“ 2 for Condition 3. In order to 
compute 6 the differential equation is used to replace the derivatives of second order in r. 
This yields: 


(36) 


*(*) = - 




+ T ±1) 


U+l) 


(37) 


b — T — 


0(z 3 ), z — * oo 

-i 
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3. THE CYLINDRICAL WAVE EQUATION: A CASE STUDY 



Fig. 1 . Errors for r — 2. Condition 1: Solid line, Condition 2: Long dash - short space, Condition 3: 
Short dash - long space , Condition 4: Equal dash and space . 


Condition 4: 


(38) 




a i \ 
Ft + Tr + s ) 



This condition is a generalization of those proposed by Barry, Bielak and MacCamy [1], 
Here 6 > 0 is a free parameter. The results are somewhat sensitive to its choice. We didn’t 
make a concerted effort to find an optimal 8, but do present the best results we obtained. 
These were for <5 = .005. Note that this condition captures the first two terms of the large 
z expansion and agrees with T at z = 0. We have: 


(39) 

*(*) = — 2 

(40) 

•-'■{S' 


z(z + \ + 8r) 
(z + 8t) 


,-U 


oo 

0 


For the numerical experiments we used second order finite difference discretizations 
with a time step of .01 and a mesh width of .02. The artificial boundary was placed at 
r = 2,4, 8. The errors were measured by comparison with a solution for t = 252 at r = 1.3 
and r — 1.8. 

The results are summarized in Figures 1-3. The maximum error for Condition 1 is 
largest, peaking at early times and slowly decaying for large time. Condition 2 is better 
for short times but has an error which grows monotonically. Condition 3 is by far the 
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best for short times but with a monotonically growing error also. Condition 4 exhibits the 
best overall performance. Its error displays oscillatory behavior. We also carried out a 
computation with mesh stretching, the results being displayed in Figure 3. We found that 
the largest mesh width had to kept fairly small to get reasonable results. Note that the long 
time errors decay slowly as r increases. The decay for Conditions 2 and 3 predicted by the 
analysis of [ 2 ] is, in particular, not evident. This is a result of the long time nonuniformity 
of the high frequency - large r asymptotics. Quantitatively, these results are predicted by 
the analysis of the behavior of 6 — T for large and small z. 

In no case do we observe a very rapid decay of the error with time. This is explained by 
the behavior of b — T for z small. As the j—j term cannot be canceled by a rational function, 
37(6 - f) will always be singular at the origin. This implies a lower bound on the decay 
in time of the inverse transform of (6 - T). The only way around this difficulty is to use 
an operator which is nonlocal in time. In joint work with Hariharan and MacCamy, the 
derivation, analysis and efficient implementation of such a boundary operator is currently 
under consideration. 
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